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The time dependent cluster approximation called the path probability method (PPM) is ap- 
plied to a pseudo-spin Ising Hamiltonian of the Slater- Takagi model for KH2P04-type hydrogen- 
bonded ferroelectrics in order to calculate the homogeneous dynamical susceptibility x{'-^) above 
and below the ferroelectric transition temperature Tc. Above the transition temperature all the 
calculations are carried out analytically in the cactus approximation of the PPM. Below the 
transition temperature the dynamical susceptibility is also calculated accurately since the ana- 
lytical solution of spontaneous polarization in the ferroelectric phase can be utilized. When the 
temperature is approached from both sides of the transition temperature, only one of relaxation 
times shows a critical slowing down and makes a main contribution to the dynamical suscepti- 
bility. The discrepancy from Slater model (ice-rule limit) is discussed in comparison with some 
experimental data. 

KEYWORDS: KDP(KH2P04), phase transition, CVM (cluster variation method), PPM (path probability method), 
dynamical susceptibility 



§1. Introduction 

Recently, we successfully applied the cluster vari- 
ation method (CVM)i) to the Slater-Takagi modeP^ 
for KH2P04(KDP)-type hydrogen-bonded ferroelectrics 
above and below the transition temperature'^-' to explain 
the anisotropy of the wave-number dependent suscep- 
tibility x(g) observed in the neutron scattering experi- 
ment.^^ On the other side, the path probability method 
(PPM)^^ devised by Kikuchi is the time dependent clus- 
ter variation method and has been applied to various 
phase transitions and transport phenomena. Its char- 
acteristic is that the stationary solution of the kinetic 
equation given by the PPM yields the equilibrium solu- 
tion obtained from the CVM in the corresponding ap- 
proximation. Further, since the PPM provides a sys- 
tematic approximation for the kinetic problem, it makes 
possible to calculate the dynamical susceptibility beyond 
the usual molecular field approximation. 

A few years ago Matsuo et al7^ re-examined the ex- 
cess entropy obtained from their own data and the other 
experimental data of heat capacity for KDP. They disc- 
cused about the discrepancy from the ice-rule of Slater 
model®) and emphasized the significance of excitation 
level in the Slater-Takagi model. 

The present purpose is to calculate the dynamical sus- 
ceptibility for KDP based on the Slater-Takagi model 
and to compare our results with experimental data of 
the dynamical susceptibility over all the temperature 
regime. Though Yoshimitu and Matubara^^ have al- 
ready calculated the dynamical susceptibility for the es- 
sentially same model for the KDP above the transition 
temperature, their calculation seems to be limited to the 
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paraelectric phase. In the present paper, not only in 
the paraelectric phase but also in the ferro-electric phase 
we calculate accurately the dynamical susceptibility for 
KDP based on the above mentioned Slater-Takagi model 
by making use of an analytical expression for the spon- 
taneous polarization. 

§2. Formulation 

There are N PO4 tetrahedra and 2N protons around 
PO4 tetrahedra in the KDP-type crystal as shown in 
Fig. 1. The pseudo-spin Ising Hamiltonian H for a con- 
figuration of 2N protons has a form^^' 



H = 



^ |^i/o(CTi, CTj, CTfe, cr/) - ^E{a^ 



(ijki) 



(2.1) 



with 

Ho{ai,aj,ak,(Ti) 



-V5{(Ji(Jk + CTjCn) - ViCTiajakCri + C 

(2.2) 

where the sum (ijkl) runs over four protons ijkl around 
each PO4 tetrahedron in the crystal, ai — ±1 stands for 
the site of the i-th proton in the double well potential 
along the 0-0 bond (hydrogen bond) between two near- 
est neighbor P04's, /id is the magnitude of an electric 
dipole moment associated with a complex K-H2PO4 and 
E is an external electric field. As is seen in Fig.l, we use 
a convention that when the i-th proton is located on the 
closer site to an O atom at the top ( bottom ) of the PO4 
tetrahedron along the easy z-axis, the i-th proton takes 
CTi = +1(— 1). The energy parameters V2,Vi,V^ and C 
are related to those of the Slater-Takagi model shown in 
Fig. 2 as 



V2 = £2/ 



Vi = -£0/4 + £1/2 -£2/ 



(2.3) 
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Fig. 1. 2-axis projection of hydrogen bonds connecting PO4 com- 
plexes and CTi , (Tj , (Tk, <Ti showing the four different pseudo-spins 

for protons around a PO4 tetrahedron. The numeral in the cen- 
ter of each PO4 tetrahedron shows relative heights of PO4 along 
z-axis. 
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Fig. 2. Energy levels, magnitude of dipole moments of a K- 
PO4 complex and probability for proton configuration in Slater- 
Takagi model. 



V5 = £0/4 - £2/8 , C = So/4: + £1/2 + £2 

where the pseudo-spin Hamiltonian (2.1) is a modified 
one used by Tokunaga et al}^^ by allowing all the con- 
figurations of four protons around a PO4 tetrahedron. 
Here the ice-rule limit characteristic of the Slater model 
for KDP is realized when (ei — £o)/£o 00 with £2 > £1. 
The ice rule is described by proton configurations in 
which (A) only one proton exists on each hydrogen bond 
between two nearest neighbor PO4 tetrahedra and (B) all 
the PO4 tetrahedra have exactly two protons adjacent to 
them. 



We now apply the path probability mcthod(PPM) in 
the cactus approximation to the present system to find 
the kinetic equation for protons. The cactus approxima- 
tion equivalent to Slater's trcatmcnt^^ takes ac(;ount of 
the proton correlations around PO4 as well as the site 
of a proton in the double well potential on each hydro- 
gen bond. However, since the derivation of the kinetic 
equation by the PPM is a little lengthy, though the final 
kinetic equation is relatively simple, here we only men- 
tion the idea of the PPM^'-'^^^. In equilibrium statistical 
mechanics, the realized state of a system in thermal con- 
tact with a heat reservoir is the minimum state of its 
free energy. When the system is not in equilibrium, we 
are interested in the time evolution of the system. The 
PPM is a method for determining the time evolution of 
the system. The idea of the PPM is to calculate a transi- 
tion probability of the ensemble of equivalent systems in 
a short time interval At from time t to t+ At. This tran- 
sition probability is called the path probability function. 
Then, we assume an extremum principle that the path 
maximizing this path probability function determines the 
time evolution of the system. 

Now, in the cactus approximation of the PPM, the 



homogeneous state of the present system at time t is 
described by five state variables defined by 

'fn{t) = {ai)t , s{t) = {(Ti<Jjak)t , qit) = {(Ti<Jj)t , 

qT>{t) = {(ri(Tk)t ) 94(i) = {(ri(Tj(Tk(Tl)^ (2.4) 

where ciach state; variable represents the correlation of 
protons ijkl around a PO4 cluster at time f(Fig.l) and 
{■ ■ \ is an thermal average at time t. After some ma- 
nipulations of the PPM we obtain a generating function 
from which a set of kinetic equations are derived through 
differentiation of interaction parameters. The generating 
function is given by 

G{L) = eTrpi{c7i,t)e-^'^'^' 

i 

jkl Pl{o-i,t) 

(2.5) 

AtHo{ai,aj,(Tk,<yi) 

= Ho(-ai, aj,ak,cri) - HQ((T.i,aj,ak,(Ti) 

where Tr^ and Tr^-fe; denote a trace operation X1<t =±i 
and J2aj,ak,cri=±i' respectively, /3 = l/fceT is the in- 
verse temperature, is a microscopic relaxation time 
of an isolated proton, Li = (5iJi^E/2 and Ai defines an 
energy increase under an inversion of only (n variable 
into —CTi. Further, pi{(Ti,t) and p4{(Ti, (Tj , (7k, <7i,t) are, 
respectively, the probability of finding the site ai of a 
proton in the i-th bond and the probability of finding 
the sites ai, crj , (Jk, cri of protons i, j,k,l around a PO4 
and are given in terms of above defined five state vari- 
ables by 

Pi{a„t) = ^{l + m{t)ai) , 
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a definition he = exp {(imE) 



P4{(Ti,crj,crk,(Ti,t) 

= ^ (l + m{t){(Ji + (Jj + (Jfe + ai) 

+q{t){aiaj + CTjCTfc + crfccr; + aicji) 
+qT){t){aiak + o-jai) 

+s{t){aiajak + crjakcri + (Tkcriai + aiaicrj) 

(2.6) 

Then a set of kinetic equations are given in a convenient 
form:^^) 



dmi{t) ^ dG{L) 
= 4 lim 



(i = 1 ^ 5) (2.7) 



dt i/3-»o dLi 

Here, it should be noted that in order to write the above 
expression an extra interaction term is virtually added 
to Hamiltonian (2.2) as 

— > Ho{ai,aj,ak,ai) 

(2.8) 

and V3 is, however, put to zero just after differentiation 
with respect to Z/3 = /3Vs in eq.(2.7). We also redefine or- 
der parameters as mi{t) = 4m{t),m2(t) = 4:q{t),m3(t) =^ 
4s(t),m4(t) — q4:{t) and m^{t) = 2q^{t) and the corre- 
sponding fields as Li{t) = 0fj,bE{t)/2,L2 = 0V2,L3 = 
pVs,L4 = pVi and L5 = /JV5, respectively. 

§3. Thermal equilibrium 

In order to obtain the dynamical susceptibility x(w) 
as the linear response to the external field, equilibrium 
values of the order parameters are required. Since the 
equilibrium state is more easily obtained from the CVM 
than from the stationary solution of the kinetic equation 
(2.7), we apply the cactus approximation of the CVM to 
the present system.^' The variational free energy F is 
obtained by 

F = U-TS (3.1) 
where the internal energy U is given by 

U/N = ^AV2q - 2V5gD - Vm - ^HAmE (3.2) 
and the entropy S is given by 
5/7V/CB =2Ti-L(ai)lnpi(ai) 



Tr 

ijkl 



P4{cri,aj,ak,crk)l^P4{(Ti,(Tj,ak,(Ti) . 

(3.3) 



The minimum condition of the variatinal free energy with 
respect to m, s, q, qn, 94 yields equilibrium relations with 



c+ _ 1 + m\ " 4 d+ _ f 1 + m \ 2 



1 — m 



d. 



1 — m 



C2 ri2 c+c_ 1 



Co ?7o co^ 770^ ' co^ \riQ 



(3.4) 



where c+, c_, d+, d_, C2 and cq shown in fig. 2 are de- 
fined from Pi{(Ji, (7j,ak, cr;) by 

1 . 

c+ = ^(1 +4m + 'iq + 2qu + 4s + 94) , 
c- = ^(1 - 4m + 4g + 2qu - 4s + 94) , 



(3.5) 



d+ = ^(1 + 2m -2s -94), 

d- = ^(l-2m + 2s-g4), 
1 . 

C2 = ^(l-4g + 2gD + ?4), 

Co = ^(1 - 2(7D + 94) 

with a normalization condition 

c++C-+4:Co + 4:{d++d-)+2c2 = l. (3.6) 

From eq.(3.4) and eq.(3.6) it is easy to solve 
c+,C-,d+,d-,C2, Co in terms of the polarization m and 
field variable h^: 



The spontaneous polarization mo^°^ is determined by 
the relation 



m = pi{+l) - Pi(-l) = c+ - c_ + 2{d+ - d-) 



as 



mo = < 



'1- 



(1 - 2770 - m? 



for T <Tc 



for T>Tc 



(3.^ 



(3.9) 



and the electric susceptibility is given by^^ 

^ _ 2(l + r7i(l-2m^)/v/r^) 

fcBT-i + 2r?o + ??2 + 2r7iMl-m^)3 ^ 

The other equilibrium order parameters in eq.(3.7) with- 
out electric field are found from eq.(3.9). Further, utiliz- 
ing eq.(3.3) the entropy versus temperature are shown in 
Fig. 3, in order to compare our result with experimental 
data of Matsuo et alJ^ In the following figures energy 
parameters such as si are measured in the units of Boltz- 
mann factor fce- 

§4. Dynamical response in paraelectric phase 

In the paraelectric phase we can carry out all the cal- 
culations analytically. In thermodynamic equilibrium of 
the paraelectric phase, under the inversion of external 
field E, the order parameters m and s are changed into 
— m and — s, respectively, while q, qr, and 54 are invari- 
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A + m Co 
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, 1 — m Co 
1 + m ryo 



?72 

C2 = Co : 

Vo 



(3.7) 



1 



m rjiCo 



1-m rjo 



d- = he 
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Fig. 3. Temperature dependence of the entropy for ei 
2£o, 6£0j lOeoi 14£o and £2 = 4£i — 2£o . 



ant. Thus m{t) and s(t) are long range order param- 
eters responding linearly to external field E{t), while 
q{t),qB{t) and q4{t) are short range order parameters 
responding quadratically to the field. Then, in order to 
obtain the linear dynamical susceptibility of the present 
system above the transition temperature, the short range 
order parameters q{t), qoit), 94 (t) can be replaced by the 
values at thermal equilibrium in paraelectric phase with- 
out electric field. A set of kinetic equations (2.7) with 
five equations is reduced to two closed equations for long 
range order parameters m{t) and s{t) up to a linear order 
to the external field by 



with r]o = e ri\= e , 772 



m(l 



e~^^2 and 



11 = 



1 + 2770 + 4771 

7?r(3- 2^- 



m 



(4.2) 



(4.3) 



1 + 2% + Arji + 772 

Now we assume fifim{t) — x{^) E eiij> icut and fids(t) = 
Xs{'-^)E cxpiiot. Substituting these relations into this set 
of kinetic equations, we finally obtain the dynamical sus- 
ceptibility x('^) by 



where A±, D are given as 



Especially the static susceptibility Xstat is obtained by 
putting a; = in eq.(4.4) as 



Xstat = X('^ = 0) 



Ml 



2(1 + m) 



(4.6) 



fefiT -1-^2770 -1- 2771 +772 

The transition temperature Tc to the ferroelectric phase 
is determined as a divergent point of the static suscepti- 
bility as 

2 g — £o/fcBTc _|_ 2 g-ei/fesTc _|_ g-£2/fcBTc — ^ (4 7) 

This expression is the one obtained by Ishibashi.-'^''^ 
When the transition temperature is approached, the first 
term of eq.(4.4) shows a critical slowing down and con- 
tributes mainly to the dynamical susceptibility. 

§5. Dynamical response in ferroelectric phase 

In the ferroelectric phase a finite spontaneous polar- 
ization nio occurs as given in eq.(3.9). When the external 
electric field E{t) is applied, the short range order param- 
eters q{t), qu(t), 94 (i) have also components proportional 
to E{t) indirectly through the spontaneous polarization 
mo. Thus we assume that 

mi{t) = 4m(f) = 4mo + xi(a;)i^expia;f 

m2{t) = 4g(f) = 4go + X2(w)£expiwi 

m3{t) = 4s(t) = 4so + X3{'^)Eexpiojt (5.1) 

7724(f) = qiit) = 94 + X4('^)-£' expiwt 

7725 (i) = 2gD(i) = + X5{'^)E expiut 

where the required dynamical susceptibility is xi^) = 

/XdXi('-^)/4 and 777,0, so, 90, (Jd and q1 arc equilibrium or- 
der parameters in the absence of external electric field. 
Then, substituting eq.(5.1) into eq.(2.7), we finally ob- 
tain a set of algebraic equations for five Xi('^) (i = 1 ^ 
5) which is read in a matrix form as 



I + M]x{uj) = b 



(5.2) 



Since the explicit forms of a matrix M and a column 
vector b are complicated and lengthy compared with the 
paraelectric phase, they are given in Appendix. Refer- 
ring to Appendix, the elements of the matrix M and the 
column vector b can be expressed in terms of only 770 , 771 
and 772 without unknown state variables since order pa- 
rameters at equilibrium are obtained analytically. Then 
the algebraic equation can be easily calculated numeri- 
cally for fixed uj using the Gaussian elimination method 
for linear algebraic equation. With the relaxation time 
Tj (1 = 1^5), the dynamical susceptibility per proton 
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dt 

ds{t) 



Li+{1- 



dt 



80Xii 



Li+ll 



2??i + 2^ + ^ -1\ 2r]i - 2^ - ^ + 1 ■ 



2r]i+2^+^-l -6?7i +2y/^+y/^+3 



+ 



2r]i - - + 1 6i]i + 2^+^ + 3 



m{t) 



(4.1) 



8^X 



s{t) 



ex fD-8A_X D-8A+X 



UbT A+ + A_ \iuj + M_A iuj + M+A 



(4.4) 



X{uj) can be written for T < Tc as 



5 



Xi 



(5.3) 



where the relaxation times are obtained by a diago- 
nalization of M in terms of a matrix U as 



{JJMU-% = ^6ii 



(5.4) 



{5ij : Kronecker's delta) 

and the intensity coefBcients Xi ^'^^ given by 

1 ^ 
Xi = ^Siidn ^(C/ ^)uUijbj . 

Especially for T > Tp. X3)X4 and Xs representing the 
intensity from each relaxation mode reduce to zero and 
in consistency to eq.(4.4) we obtain 



Xi 



+ ioJTi 



(5.5) 



Though there appear five relaxation times in the fer- 
roelectric phase, only one of them shows the critical 
slowing down and contributes mainly to the dynam- 
ical susceptibility when the transition temperature is 
approached (Fig. 4). The real X {^) a-nd imaginary 
X (w) part of the dynamical susceptibility are defined as 
X{^^) = x'(w) - ix"{^^)- The results for x'(w) and x"(w) 
versus to and temperature (T — Tc)/Tc in the para- and 
ferro-electric phases are shown in Fig. 5. 

§6. Results and discussions 

The result shows that, though there are, respectively, 
two and five relaxation times in para- and ferro-electric 
phase according to a set of independent kinetic equa- 
tions, only one of the relaxation times shows a criti- 
cal slowing down when the temperature T approaches 
the transition temperature Tc from above and below the 
transition temperature and makes a main contribution 
to the dynamical susceptibility. 

In order to compare the experimental data with our 
results we present the; temperature dependence of the 
real part x'i'^) imaginary part x"{'^) ^ov constant w 
for various values si (Fig. 6). The hilly behaviors appear 



61=813 K , rc=123 K 
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o 



tr 




-0.4 -0.2 
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T-T 
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Fig. 4 
T > 



Relaxation times versus temperature for ei = 813 K. For 
Tc, only Ti and T2 take part in x{<^)- 



not only above the transition temperature Tc but also 
below Tc for the finite ei whereas in the ice-rule limit 
X(w) vanishes completely below Tc. The dip of x (f^) at 
T = Tc is caused by the vanishing of the numerator due 
to the contribution from the relaxation mode showing a 
critical slowing down. The experimental data^^^ show 
the hilly behavior below the transition temperature and 
the dip structure at the transition temperature. In our 
calculation the contribution from the relaxation mode 
showing a critical slowing down overwhelms contribu- 
tions from other modes and the dip goes to almost zero 
contrary to the experimental data. Recently. Matsuo et 
al7^ re-examined experiments of the heat capacity and 
estimated the transition entropy AS due to proton order- 
ing from the experimental data. They discussed the dis- 
crepancy from the Slater theory and estimated the con- 
tribution from the excitation level of the Slater- Takagi 
model. We presented the entropy curve versus tempera- 
ture from our calculation for various parameters in fig. 3. 
These results reveal that the ice-rule in the Slater model 
is not completely satisfied in the KDP crystal. 

The dynamical susceptibility has been calculated not 
only above the transition temperature but also below it 
in the cactus approximation of the Slater- Takagi model 
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Fig. 5. w and (T - rc)/Tc dependence of x'ii^) and x"('^) for 
£0 = 85.6K, £1 = 813K, £2 = 4£i - 2£o. 



utilizing an analytical solution for the spontaneous po- 
larization. The results based on the Slater- Takagi model 
are in good agreement with the experiments on dynam- 
ical susceptibility and excess entropy. 
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Appendix: Derivation of eq. (5.2) 

The generating function (2.5) is conveniently written 



as 



G(Z) = 0Tr 



1 + m{t)(Ji \jki 
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Fig. 6. Temperature dependence of the susceptibility x('^) foi' 
various £1 under the constant w. The solid and dotted curve de- 
notes, respectively, the real x'i'^) and the imaginary part x"{<^) 
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where L is the energy parameter vector defined as 

L = p{C V2 Vs V4 V5) 

and the Ising spin vector a is defined as 
/ 1 

■ ■ o-j + o-fe + (Jl 

iCFj + O-jCTfe + (JfcCr; + aiCFi 



(A-2) 



(Tfe + (TjCTl 



(A-3) 



The thermal average of the Ising spin vector is the order 
parameter vector defined by 

m{t) = (1 4m(t) Aq{t) 4s(t) q^t) 2qr>{t)) . 

(A-4) 

The probability P4:{{<y}ijki){= Pii^^i, o-j, cfk, o^i)) can 
also be written by using cr and m(t) as eq.(2.6). Then 
a set of kinetic equation (2.7) is written in a vector form 

by 

dAm^) _ dG{L) 

and is further rewritten to a linear order of the external 
electric field E as 



Edited by J.L. Moran-Lopez and J.M. Sanchez, Plenum Press, 
New York, 1996. 

13) K.Wada and M.Kaburagi, Prog.Theor.Phys.Suppl. 115 293 
(1994). 

14) G.V.Kozlov, S.P.Lebedev, A.A.Minaev, A.A.Volkov and 
E.V.Monia, Ferroelectrics 2 1 373 (1978). 



where p'i{a) = (1 + mQa)/2 and P%{{(y)ijki) = 
p1{(Ti,(Tj,ak,(Ti) are thermal equilibrium probabilities 
without external field E = and Am{t) is a linear dif- 
ference from equilibrium value of order parameter vector 
m{t). The variable h and the vector h' are further de- 
fined as 

^{{<y}i3u) = h'{ui,uj,ak,cji) = e^'^-^^'^'^Aia- 

where Lq is an interaction parameter vector L with E = 

0. In order to obtain the dynamical susceptibility x(cj), 
Am{t) = x{^)E expiijot is assumed and then eq.(A-6) 
can be easily rewritten into a form (5.2) in the main 
text: 

(^7 + M)xH = b- (A-7) 
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